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ABSTRACT 

We analyze a sample of optical light curves for 100 quasars, 70 of which have black hole mass 
estimates. Our sample is the largest and broadest used yet for modeling quasar variability. The 
sources in our sample have z < 2.8, lO'^^ < AL;,(5100A) < 10'*^, and 10^ < Mbh/Mq < 10^°. We 
model the light curves as a continuous time stochastic process, providing a natural means of estimating 
the characteristic time scale and amplitude of quasar variations. We employ a Bayesian approach to 
estimate the characteristic time scale and amplitude of flux variations; our approach is not affected by 
biases introduced from discrete sampling effects. We find that the characteristic time scales stongly 
correlate with black hole mass and luminosity, and are consistent with disk orbital or thermal time 
scales. In addition, the amplitude of short time scale variations is significantly anti-correlated with 
black hole mass and luminosity. We interpret the optical flux fluctuations as resulting from thermal 
fluctuations that are driven by an underlying stochastic process, such as a turbulent magnetic field. 
In addition, the intranight variations in optical flux implied by our empirical model are < 0.02 mag, 
consistent with current microvariability observations of radio-quiet quasars. Our stochastic model 
is therefore able to unify both long and short time scale optical variations in radio-quiet quasars as 
resulting from the same underlying process, while radio-loud quasars have an additional variability 
component that operates on time scales < 1 day. 

Subject headings: accretion, accretion disks — galaxies: active — methods: data analysis — quasars: 
general 



1. INTRODUCTION 

It is widely accepted that the extraordinary activity 
associated with quasars^ involves accretion onto a super- 
massive black hole, with the UV/optical emission arising 
from a geometrically thin, optically thick cold accretion 
disk. Aperiodic variability across all wavebands is ubiq- 
uitous in AGN, with the most rapid variations oc curing 
in the X-rays (for a review, see lUhich et aLlHool . The 
source of quasar variability is unclear, and several models 
have been proposed for describing the optical variabil- 
ity of quasars, includin g accretion disk i nstabilities (e.g., 
Kawaguchi et al.iri998D. supe rnovae fe.g.. lAretxaga et all 
1997( ). microlensing |llawki ns 2000), and more genera l 
Poisson process models (e.g., Cid Fcrnan des et al.ll2000f ). 
However, recent results from reverberation mapping have 
shown that the broad emission lines respond to varia- 
tions in the continuu m emission after some time lag (e.g., 
iPeterson et al.ll2004f ). implying that the continuum vari- 
ations are dominated by processes intrinsic to the accre- 
tion disk. If the optical/UV variations are intrinsic to the 
accretion disk, then thermal fluctuations appear to be a 
natural choice for driving the optical/UV variations, as 
the optical/UV emission is thought to be thermal emis- 
sion from the accretion disk. The fact that quasars be- 
come bluer as they brighten is con sistent with a ther- 
mal origin (e.g., IGiveon et all 119991 : iTrevese et all 120011 : 
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iGeha et al.ll2003 h. Moreover, the thermal timescale is 
sensitive to the disk viscosity. Variability is therefore a 
potentially important and powerful probe of the quasar 
central engine and accretion disk physics. 

A considerable amount of our interpretation and un- 
derstanding of AGN accretion disks, and consequently 
their optical/ UV emission, is based on t he so-called a- 
prescription (jShakura fc SyunvaevI Il973h . Within the 
standard a model, the viscosity, thought to be the source 
of the thermal disk emission and outward transfer of 
angular momentum, is parameterized as being propor- 
tional to the total pressure in the disk. Previous work 
on estimating a from quasar variability has found a 
value of g ^ 0.01 (e.g., ISiemiginowska fc Czernvl 119891 : 
ICollier fc Peterscml I2OOII : IStarling et all l200l . How- 
ever, when the disk is dominated by radiation pres- 
sure, as is thought to be the case in the inner re- 
gions, an g-disk is bo th ther mally and viscously u nstable 
(jShakura fc SunvaevI [l97l iLightman fc EardlevI [l97l . 
In particular, for a radiation pressure dominated disk, 
the thermal instability is expected to grow exponentially 
on a time scale similar to the thermal time scale. For 
AGN, the thermal time scale is on the order of months 
to years, and in general there is no evidence for insta- 
bilities in t he optical light curves o f AGN which span 
> tth (e.g.. IColher fc Peterson|[2001h . However, the op- 
tical light curves for a few sources may s how evidence 
of instability on a thermal ti me scale (e.g.. ICzernv et aTl 
[200llLub fc de Rmtejll99l l. 

If the stress is not proportional to the total pres- 
sure, but rather some other combination of the gas and 
radiation pressure, then the dis k becomes more stable 
to thermal p erturba tions (e.g.. IStella fc Rosnerj 119841: 
Szuszkiewic? 119901: iMerloni fc FabianI 120021 : iMerlonil 
2003| ). Alternatively, if part of the accretion energy 
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is dissipated in a hot corona, then there is no longer 
any runaway heating in the dis k, as additional cool- 
ing i s provided by the corona (jSvensson fc Zdziarskil 
Il994f) . For example, disk models with alternative pre- 
scriptions for a, as well as the addition of a corona 
and jet, have been able to reproduc e the variability 
of the microquasa r GR S 191 5-H05 (|Navakshin et all 
[2000t iJaniuk et alJ [2000l . [200^ . Previous work on 3- 
dimensional magneto-hydrodynamic (MHD) simulations 
of radiation pressure-dominated AGN accretio n disks 
have not observed a th ermal instability (e.g., iTurneil 
l2004t iHirose et alll2008| ). Moreover, the most promis- 
ing physical mechanism behind the vi scous torque is the 
magneto-ro tational instability (MRI, iBalbus fc HawlevI 
llQQll Il998f ). and recent numerical and analytical work 
has suggested that the a prescription m ay be a poor 
repre sentation for MRI-driven viscosity (jPessah et al.l 
|2008[) . Variability studies can therefore lend observa- 
tional insight into the inadequacies of the traditional a- 
prescription, and possibly lead to a more accurate charac- 
terisation of the relationship between stress and pressure 
in the accretion disk. 

A successful model for quasar X-ray variability de- 
scribes the X-ray variations on long time scales as being 
the result of perturbations in the accretion rate that oc- 
cur o u tside of the X-ray emitting region (e.g., Ly_ub arskii 
[19971 : iMaver fc PriMy 120061 : iJaniiik fc CzeT^ 120071 ). 
These accretion rate perturbations then travel inward, 
modulating the X-ray emitting region. It has been sug- 
gested that the origin of such perturbations is the result 
of a magnetic field randomly varying in tim e, and may 
be re lated to the appearance of an outflow ([King et all 
|2004[ ). If this model for the X-ray variability is correct, 
we would expect to also see variations in the optical lumi- 
nosity, whose origin lies in the disk at radii farther from 
the central source. Therefore, understanding the origin 
of quasar optical variations will not only lead to a bet- 
ter understanding of accretion disk physics, but may also 
lead to a better understanding of the origin of quasar X- 
ray variability, possibly unifying the source of variability 
in the two bands. 

Quasars have also been observed to vary on time 
scales as short as hours (so-c alled 'microvariabil- 
ity' o r 'intranigh t var i ability ', iGopal-Krishna et all 
2001 iStalin et aLl [200l [20051: iGupta fc Joshil 120051 : 
Carini et al.l 1200^ . Microvariability is known to be 
stronger in radio-lou d quasars, especially blazars (e.g., 
iGupta fc Joshil l2005f l. while microvariability in radio- 
quiet quasars is usuall y not detected above the photomet- 
ric u ncertainty (e.g.. IGupta fc Joshil l2005t ICarini et al.l 
l2007t l. The enhanced level of microvariability in radio 
loud objects suggests that it is due to processes in a jet, 
while the physical cause of microvariability in radio quiet 
objects has remained a puzzle. Reprocessing of X-rays, 
disk instabilit ies, or a weak bl azar component have been 
suggested (Cz ernv et alll2008[) . 

There have been numerous previous investigations of 
quasar optical variability. However, because of the diffi- 
culty in obtaining high quality, well sampled light curves 
that cover a long time span, most previous work has in- 
volved ensemble studies of quasars, or analysis of sim- 
ple correlations involving variability amplitude. The 
most well known result from previous work is a ten- 
dency for AGN to become less variable as their luminos- 
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have also b een claims of a vari a bility-redshift corre- 
lation (e.g., ICr istiani et al. 1990'; 'Cid Fernande s et al.l 
119961 : ICristiai^et al. 1996: Trevese fc Vagnetti.200l l. al- 
though the sign of this c orrelation varies bet ween stud- 
ies (e.g., see the list in iGiveon et al.l Il999l ). If real, 
the variability-redshift correlation is most likely caused 
by the fact that quasars are more variable at shorter 
wave l engths (e.g.. ICutri et al.lll985l: Idi Clemente et al.l 
[l99^lHelfand et al.ll2001l : IVanden Berk et al.ll2004D . cor- 
responding to regions in the disk closer to the central 
black hole. Recently, a correlation between optical vari- 
abilit y and black hole mass has been claimed (jWold et al.l 
|2003). While many of these correlations are formally sta- 
tistically significant, they often exhibit considerable scat- 
ter when one measures variability for individual objects. 

A few previous studies have employed spectral tech- 
niques, such as power spectra and structure func- 
tions, in the analysis of quasar optical light curves. 
From these studies it has been inferred that quasar 
optical light curves generally have variations of 
10% on timescales of months, and that the power 
spectra of optica l li ght curves_ is well described as 
P(f) PC 1/f^ (| Giveon et a l. i9_99:; Collier fc Peterso3 
1200 ll : ICzernv et al.ll2063h ; power spectra of this form are 
consistent with random walk, or more generally, autore - 
gressive processes. In addition [Collier fc PetersonI ()2001l ) 
analyzed a sample of optical light curves from 8 low-z 
Seyfert 1 galaxies. They found that the characteristic 
time scales of optical variations for the AGN in their 
sample are ~ 10-100 days and correlate with black hole 
mass, consist e nt wit h disk orbital or thermal time scales. 
ICzernv et al.l ()1999f ) found evidence for a flattening of 
the optical power spectrum of NGC 5548 on time scales 
longer than ^100 d ays. Both ICzernv et all (i_1999) and 
ICzernv et al.l ()2003[ ) also suggested that the long time 
scale optical variability was due to thermal fluctuations 
in the accretion disks of NGC 5548 and NGC 4151, re- 
spectively. 

Motivated by the potential in variability studies for in- 
creasing our understanding of the structure of quasar ac- 
cretion disks, we have compiled a sample of well-sampled 
optical light curves from the literature. We directly 
model the quasar optical light curves as a stochastic pro- 
cess, in contrast to previous work based on more tradi- 
tional Fourier (i.e., spectral) techniques. Our method al- 
lows us to describe quasar light curves with three free 
parameters: a characteristic time scale, amplitude of 
short time scale variability, and the mean value of the 
light curve. In addition, our method enables us to es- 
timate the characteristic time scale of quasar variations 
without the windowing effects that can bias spectral ap- 
proaches. Our sample consists of 100 AGN, 70 of which 
have black hole mass estimates. The sources in our 
sample have z < 2.8, 10^^ < AiA(5100A) < 10^^, and 
10^ < Mbh/Mq < 10^°, making this by far the largest 
sample yet used for this kind of study. 

In this work we adopt a cosmology based on the 
WMAP results (h = 0.71, = 0.27, f^A = 0.73, 
ISpergel et an[2003l) . 

2. DATA 
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In this work we analyze a sample of 100 quasar op- 
tical light curves, compiled from the literature. Our 
sample consists o f 55 AGN from the MACHO survey 
(|Geha et al.l l 2003l). 37 Palom ar Green (PG) quasars from 
the sample of iGiveon et alJ (1999), and 8 Seyfert galax- 
ies from the AGN Watch^ database. We were able to 
obtain black hole mass estimates for 71 of the AGN, 
where Mbh has been estimated for 20 of th em from re- 
verberation mapping (jPeterson et al.|[2004D . and Mbh 
is estimated for the remaining 51 AGN from the broad 
emission lines using standard scaling relationships (e.g., 
IVestergaard fc Petersonll2006f l. Our sample is summa- 
rized in Table [H 

2.1. Macho Quasars from Geha et al.(2003) 

We collected R band light curves from AGN se- 
lected by the MACHO su rvey based on their variability 
(|Alcock et all 119971 [T999I) . The motivation for the MA- 
CHO survey was to study Galactic microlensing events 
behind the Magellanic Clouds. However, the survey was 
also able to select quasars via their variability, confirmed 
via spectroscopic follow-up, producing 59 quasars with 
well-sampled light curves over abroad range in redshift 
(0.1 < z < 2.8) (iGeha et al.ll2003f ). The quasar light 
curves span ^ 7.5 years, and have a typical sampling 
interval of ^ 2-10 days, although much longer gaps ex- 
ist for some light curves. In general, the MACHO light 
curves are the highest quality in our sample, often be- 
ing frequently and regularly sampled. However, we note 
that because the MACHO sources were selected based 
on their variability, this sample is bia sed toward more 
variable objects. See iGeha et al.l (|2003D for more details 
of the MACHO quasar catalogue. 

Spectra for the MACHO quasars are presented in 
I Geha et al] ()2003[ ). and spectra for 27 of these quasars 
were kindly provided to us by Maria Geha. We calcu- 
lated black hole mass estimates from the source luminos- 
ity and the FWHM of the H^, Mg II, or C IV broad 
emission line using stan dard scaling relationships (e.g., 
IVestergaard fc Peterson 2006, Vestergaard 2008). How- 
ever, the spectra are not flux-calibrated, so the luminosi- 
ties at 1350, 3000, and 5100A were estimated from the 
photometric data, assuming a power-law continuum with 
a = 0.5,/^ oc (jRic hards et al. 2001). Typical uncer- 
tainti es on the broad line mass estim ates are ~ 0.4 dex 
(e.g.. IWstergaard fc PetersonI l2006l ) . When combined 
with the measurement error in the FWHAI measure- 
ments, and the uncertainty due to the intrinsic scatter in 
quasar spectral slopes (~ 0.3, Richards et al. 2001), our 
typical adopted uncertainty in Mbh was ~ 0.45 dex for 
the MACHO quasars. 

2.2. PC Quasars from Civeon et al.(1999) 

iGiveon et al.l (|1999f) report B and R photometric light 
curves for 42 PG quasars spanning ~ 7 years with a 
typical sampling interval of ^ 40 days. These quasars 
are all bright ( i3 < 16 mag) and nearby (z < 0.4). 
iPeterson et al.l (|2004D report estimates of Mbh and 
La(5100A) calculated from reverberation mapping for 12 
of these qua sars. Estimates of Mbh and La (5100A) were 
taken from IVestergaard fc PetersonI (120061 ) for the re- 

^ http:/ /www. astronomy.ohio-statc.edu/~agnwatch/ 



mammg IGiveon et al] (fT999l) quasars. See IGiveon et all 
(|1999D for further details. 

2.3. Seyfert Calaxies from ACN Watch Database 

The light curves for the remaining 8 AGN in our sample 
are from the AGN watch project. The optical light curves 
for the Seyfert Galaxies AKN 564, FairaU 9, MRK 279, 
MRK 509, NGC 3783, NGC 4051, NGC 4151, NGC 5548, 
and NGC 7469 were taken from the AGN watch website. 
In general, we used the light curves at 5100A, with the 
exception of AKN 564, for which we used the i?-band 
light curve. We excluded 3C390.3 from our analysis be- 
cause it is a classified as an optically-violent variable, and 
thus may experience a different variability mechanism 
than the other AGN in our sample. Bla ck hole masses 
and 5 IOOA luminosities were taken from IPeterson et al.l 
(|2004f l. Further details may be found in the references 
listed in Table [H Although the AGN Watch sample is 
small compared to the other two, these sources are par- 
ticularly important in our analysis because they help to 
anchor the low-L and \ow-Mbh end of the correlations 
analyzed in § [H 

3. THE STATISTICAL MODEL: CONTINUOUS 
AUTOREGRESSIVE PROCESS 

The previous analysis of IGiveon et al.l ()1999t ) and 
ICollier fc PetersonI (|2001[ ) suggests that the optical vari- 
ability can be well-described by a power-law power spec- 
trum with slope ~ 2. The power spectra of the MACHO 
quasars have not been analyzed for individual objects, al- 
though Hawkins (2007) investigated the power spectra of 
the ensemble of objects. In Figure[T]we show the geomet- 
ric mean i?-band power spectra for the MACHO quasars, 
along with the 90% confidence region on the geometric 
mean. The power spectra for the MACHO sources are 
well described by a power law of slope ~ 2, consistent 
with previous work. The lack of any peaks in the power 
spectra, as well as the aperiodic and noisy appearance of 
quasar light curves, suggests that quasar light curves are 
stochastic or chaotic in nature. 

In this work we model quasar light curves as a con- 
tinuous time first order autoregressive process (CAR(l)). 
Power spectra of the form P{f) cx 1//^ are consistent 
with a first order autoregressive (AR(1)) process. We 
model this stochastic process in continuous time, both 
because the actual physical processes in the accretion 
disk are continuous, and because doing so allows a nat- 
ural way of handling the irregular sampling of our light 
curves. Moreover, this process is well-studied, only has 
three parameters, and provides a natural and consistent 
way of estimating a characteristic time scale and variance 
of quasar light curves. 

The CAR(l) process is describe d by the following 
stoch astic differential equation^ (e.g.. lBrockwell fc Davis! 
I2OOI : 

dX[t) = --X{t)dt + <7\/lite{t) +h dt, T,a,t>0. (1) 

T 

Here, r is called the 'relaxation time' of the process X{t), 
and e{t) is a white noise process with zero mean and 

Strictly speaking, the stochastic differential equation is com- 
plicated by the fact that white noise does not exist as a derivative 
in the usual sense. However, we ignore the mathematical techni- 
calities for ease of interpretation of Equation JTJ 
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TABLE 1 
Quasars Analyzed in this Work 
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References. — (1) Giveon et al. 1999 (2) Geha et al. 2003 (3) Stirpe et al. 1994 (4) 
Peterson et al. 2000 (5) Kaspi et al. 1996 (6) Santos-Lleo et al. 2001 (7) Peterson et al. 2002 
(8) Carone et al. 1997 (9) Shemmer et al. 2001 (10) Collier et al. 1998 

Note. — Table [1] is published in its entirety in the electronic edition of the Astrophysical 
Journal. A portion is shown here for guidance regarding its form and content. 
^ Icr uncertainty on log Mbh-^ Reference for the optical light curve data. 
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1. — The geometric mean power spectrum (solid noisy line) 
for the i?-band light curves of the MACHO quasars, along with 
90% confidence region (shaded region). The thick solid line is a 
power spectrum of the form P{f) oc 1//'^ with an additive measure- 
ment error contribution. The optical light curves for the MACHO 
quasars are well described by a 1//^ power spectrum, consistent 
with other samples of quasars. Power spectra of the form 1//^ are 
suggestive of random walk and related stochastic processes. 



variance equal to one. Within the context of this work, 
X{t) is the quasar flux. We assume that the white noise 
process is also Gaussian. The mean value of X{t) is br 
and the variance is Further details on the CAR(l) 

process are described in the Appendix. 

The relaxation time, r, can be interpreted as the time 
required for the time series to become roughly uncorre- 
lated, and a can be interpreted as describing the vari- 
ability of the time series on time scales short compared 
to r. Within the context of this work, X{t) is the quasar 
light curve. It is tempting to associate r with a charac- 
teristic time scale, such as the time required for diffusion 
to smooth out local accretion rate perturbations, and a 
to represent the variability resulting from local random 
deviations in the accretion disk structure, such as caused 
by turbulence and other random magneto-hydrodynamic 



l + (27rT/)2- 

From Equation ([2]) we infer that there are two impor- 
tant regimes for Px{f): Px{f) oc 1/p for / > (2^t)-i 
and Px{f) oc constant for / < (27rr)~^. Therefore, 
the CAR(l) process has a power spectrum that falls off 
as 1//^ at time scales short compared to the relaxation 
time, and flattens to white noise at time scales long com- 
pared to the relaxation time. Because 'characteristic' 
time scales of quasar light curves are often deflned by a 
break in the power spectrum, this is an additional justi- 
fication of associating r with a characteristic time scale. 
In addition, because the power spectra of quasar opti- 
cal light curves are well described by Px{f) oc 1//^, it 
suggests that a CAR(l) process should provide a good 
description of the light curves, with r being on the order 
of the length of the light curves or longer. 

To illustrate the CAR(l) process, we simulate four 
CAR(l) light curves. The light curves were simulated 
by first simulating a random variable from a normal dis- 
tribution with mean rb and variance tct^/2; note that 
this is the mean and variance of the CAR(l) process. 
Then, from this random initial value, we simulated the 
rest of the light curve using Equations (|A4[) and (|A5|) 
in the Appendix. These simulated light curves span a 
length of 7 yrs and are sampled every 5 days. The sim- 
ulated light curves span a period in time similar to the 
quasar light curves analyzed in this work, but are bet- 
ter sampled than most of the quasar light curves. Three 
characteristic time scales of interest for quasars are the 
light crossing time, the gas orbital time scale, and the 
accretion disk thermal time scale. These time scales are 
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where Mbh is the mass of the black hole, R is the 
emission distance from the central black hole, Rs = 
2GMbh/c^ is the Schwarzschild radius, and a is the 
standard disk viscosity parameter. For the simulated 
quasar hght curves, we use Mbh — lO*M0,a = 0.01, 
and R = lGGi?s, and set r equal to each of these three 
time scales. In addition, we use b = and a = 1. Note 
that the assumed value of a only affects the thermal time 
scale, and a higher value of a results in a shorter time 
scale. 

The simulated light curves are shown in Figure [U and 
their corresponding power spectra are shown in Figure 
[21 The increased amount of variation on long time scales 
with increasing r is apparent. In addition, because ti^ 
is shorter than the time sampling, the first simulated 
light curve is only sampling frequencies on the flat part 
of the power spectrum, giving it the appearance of white 
noise. In contrast, the two simulated light curves with 
the longest time scales are sampled on the 1/ P part of 
the power spectrum, giving them more of a 'red noise' 
appearance. In addition, 'red noise' leak affects the es- 
timated power spectrum of the light curve with r — 4.6 
yrs, evidenced by the constant offset between the true 
power spectrum and the estimated one. Red noise leak 
occurs when power from time scales longer than the span 
of the time series 'leaks' into the shorter time scales, bi- 
asing the power spectrum when estimated as the mod- 
ulus of the discrete Fourier transform (e.g.. Ivan der KlisI 
fl99l . 

3.1. Estimating the Parameters of a CAR(l) Process 

The parameters for a CAR(l) process are commonly 
estimated by maximum-likelihood directly from the ob- 
served time series. This is an advantage over non- 
parameteric approaches, such as the discrete power spec- 
trum or the structure function. The observed power 
spectrum and structure function can both suffer from 
windowing effects caused by the finite duration and sam- 
pling of the light curve, whereby power from high fre- 
quencies can leak to low frequencies (aliasing), and power 
at low frequencies can leak to high frequencies (e.g., red 
noise leak). For ground based optical observations, an 
additional complication is the periodicity enforced in the 
sampling caused by the Earth's rotation around the sun, 
as objects are only observable during certain times of the 
year. For example, this periodic sampling can be seen in 
the light curve for the MACHO source shown in Fig- 
ure m All of these effects can bias the power spectrum 
or structure function when estimated directly from the 
light curve in a non-parameteric fashion. Similarly, these 
effects can bias parameter estimates when fitting a model 
to the observed power spectrum or structure function. 

In contrast, estimating a 'characteristic' time scale and 
variance directly from the observed time series, instead 
of from the observed power spectrum or structure func- 
tion, has the advantage of being free of windowing effects, 
giving unbiased estimates of r and cr^. Of course, this 
requires one to assume a parameteric model for the time 
series (or power spectrum), but as we will show below 
the CAR(l) process provides a good description of most 
of the AGN light curves analyzed in this work. Further- 




FlG. 2. — Light curves simulated from a CAR(l) process for three 
different characteristic time scales, assuming typical parameters 
for quasars {Mbh = 10^ Mq,Rs = 100, a = 0.01, see Eq.[3]-[5]). 
From top to bottom, these are the light crossing time, t = 1.1 days, 
the disk orbital time scale, r = 104 days, and the disk thermal time 
scale, T = 4.6 yrs. The stochastic nature of the CAR{1) process is 
apparent, and the light curve exhibits more variability on longer 
time scales as the characteristic time scale increases. 

more, higher-order terms can be add ed to Equation ([U) 
to al low additional flexibility (e.g., iBrockwell fc Davia 
l2002f ). but this is beyond the scope of the current work. 

When the light curve is measured with error, 
the likelihood function can be calculated using a 
'state-space' re presentation of the time series (e.g., 
IBrockwell fc Da vis 2002) . Denoting the measured fluxes 
as xi, . . . , Xn, observed at times ii, . . . , i„ with measure- 
ment error variances ct^, . . . , cr^, the likelihood function, 
p{xi, . . . ,Xn\b, cr, r), is a product of Gaussian functions: 



p{xi 



\b,a,T)^l[[27T{n, + af 

1 (Xi 



'1/2 



1=1 

X exp 
bT 



2 n, 



xo = 



(6) 

(7) 
(8) 

(9) 
(10) 

(11) 
(12) 

The maximum-likelihood estimate is then found by max- 
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Fig. 3. — Power spectra for the simulated CAR,(1) light curves 
shown in Figure [2] The actual power spectra are shown with a 
solid line, and the empirical power spectra estimated directly from 
the light curves are the noisy curves. The power spectra are flat 
on the 'white noise' part of the curve, corresponding to frequencies 
/ < {27tt)~^, and fall off as 1//^ on the 'red noise' part of the 
curve, / > {2ttt)~^. As t increases, the break in the power spectra, 
marked with a vertical line, shifts toward smaller frequencies. For 
the CAR(l) process with r = ttf^, red noise leak biases the power 
spectrum estimated directly from the simulated light curve. 

imizing Equation ^ with respect to b, t, and a. 

In this work we employ a Bayesian approach in or- 
der to directly compute the probability distribution of 
b, r, and cr, given our observed light curves. The proba- 
bility distribution of the parameters, given the observed 
data (i.e., the posterior distribution), is calculated as the 
product of the likelihood function with a prior probabil- 
ity distribution. In this work we assume a uniform prior 
on b and a. For deriving a prior on r, we note that when 
the data are regularly sampled, the CAR(l) process re- 
duces to the AR(1) process described by Equation (jAl[) 
in the Appendix with aAR = e"^*/"^, where At is the 
time sampling interval. In this work we assume a uni- 
form prior on aAR- For an AR(1) process, oar gives 
the correlation between Xi and Xi-i. Therefore, we con- 
sider it reasonable to assume that any value of aAR is 
a priori likely, i.e., we do not assume anything a priori 
about the correlations between subsequent data points, 
and therefore we assume a uniform prior on aAR from 
to 1. This prior is non-informative in the sense that all 
of the information on aAR comes from the data. 

The accuracy of the fit can be assessed by comparing 
the residuals of the light curve with the values expected 
under the assumption of a CAR(l) process. Equation 
^ implies that if the CAR(l) process provides a good 
model of the observed data, then the residuals should be 
uncorrelated and follow a normal distribution: 

^Z^ ^jV(0,l). (13) 



Here, the notation x ^ A^(0, 1) means that x is dis- 
tributed according to a normal distribution with mean 
equal to zero and variance equal to one. The goodness 
of fit can then be assessed by inspecting a plot of the 
residuals with time to ensure that they arc uncorrelated, 
and by comparing a histogram of the residuals with the 
expected standard normal distribution. 

3.2. Fitting the Quasar Light Curves 

In this work we model the logarithm of the flux as 
following a CAR(l) process, or equivalently the apparent 
magnitudes. We do this because the assumption of a 
Gaussian white noise process in Equation ([T]) produces 
both positive and negative values of X{t), while flux is a 
strictly positive quantity. The logarithm function maps 
a strictly positive quantity to the interval (— oo, oo), and 
therefore Equation ([1]) is likely to be a better description 
of the light curve for the source apparen t magnitude, as 
opposed to the source flux. Uttle v et al.l |2005) have also 
argued for describing accreting black holes as a Gaussian 
stochastic process in the logarithm of the flux. 

An additional correction is needed to correct for cos- 
mological time dilation. Because the quasars are fit in 
the observed frame, and time scales decrease as (1 + z), 
spurious correlations may arise if one does not correct to 
the quasar rest frame. This is particularly problematic 
when dealing with flux limited samples, which create an 
artificial correlation between z and both luminosity and 
Mbh- Noting that dtobs = (1 + z)dtrest, Equation ([T]) 
can be expressed in the forms 

dX{t) = —X{t)dtobs + Oobs \/ dtobs^{t) + bobs dti0) 

^obs 

= X{t)dtrest + (Jobs V(l + Z)dtrest^{t) 

^obs 

+ z)bobs dtrest (15) 

— X {t)dtrest + Crrest\'^'dtrest^{t) 

^rest 

-^brest dtrest- (16) 

From Equations (|14 p - (fT6|) it is apparent that the ob- 
served and rest frame parameters are related as 

Trest^i^ + zy^Tobs (17) 
arest = {l+z)^^^aobs (18) 
brest^i^ + z)bobs- (19) 

Noting that the mean of the CAR(l) process is br, and 
that the variance is Tcr^/2, Equations pT|) - P^ imply 
that the mean and variance of a CAR(l) process are 
unaffected by cosmological time dilation. However, the 
variance observed from a light curve with a finite dura- 
tion and sampling is still affected by time dilation, as the 
observed variance over a time interval At is the integral 
of Equation ([2]) over that time interval. In what follows, 
the quantities t, a, and b will always refer to the quasar 
rest frame quantities, unless specified otherwise. 

Random draws of b, r, and a from the posterior prob- 
ability distribution are obtained using a Metropolis- 
Hastings algorithm. We calculated an estimate of each 
parameter as the median of the posterior, and the pos- 
terior medians, standard deviations, and confidence in- 
tervals are computed using these random draws. In gen- 
eral the posterior median values were not significantly 
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Fig. 4.— Light curves and best fit CAR(l) processes for NGC 5548 (left), MACHO source 013.05717.0178 (middle), and PG 0804+761 
(right). The top panels show the light curves (data points) along with the best fit CAR(l) process (solid line), the middle panels show 
the standardized residuals (Eg. 1131 ) as a function of time (data points) and a moving average estimate (dashed line), and the bottom 
panels compare a histogram of the standardized residuals with the expected standard normal distribution. The AGN light curves are well 
described by a CAR(l) process with characteristic time scales t = 214, 6026, and 1148 days for NGC 5548, MACHO 013.05717.0178, and 
PG 0804+761, respectively. 



different from a maximum likelihood fit. The goodness 
of fit for the light curves was determined by examining 
a histogram of the residuals and a plot of the residuals 
against time, as described in § 13.11 Occasionally out- 
lying values of the flux are present for the light curves 
with more data points, possibly due to unidentified sys- 
tematic error. These outlying data points were removed 
and the light curves were refit. In Figure 2] we show 
the light curve and best fit CAR(l) model for the most 
densely sampled object in our sample, NGC 5548, for a 
representative light curve from the MA CHO sample, and 
for a repre sentative light curve from the PG sample of 
iGiveon eTal . (1999). In general, the CAR(l) model pro- 
vided a good fit to the quasar light curves analyzed in 
this work, and we only flagged 9 out of 109 light curves as 
having a bad fit. The best fit CAR(l) parameters, along 
with their uncertainties, for the 100 AGN in our sample 
are listed in Table [H The 9 objects for which the fit 
was deemed unacceptable are not used in the regression 
analysis. 

4. RESULTS 

The distribution of cr, r, and light curve standard de- 
viations for our sample are shown in Figure El The light 
curve standard deviation is calculated as the square root 
of the light curve variance, s = u^t 12. The best fit 
quasar relaxation times have a median value of 540 days 
and a dispersion of 0.64 dex, and show typical long time 
scale optical variations of 3-30 per cent. However, the 
uncertainties on r are large and make a considerable con- 
tribution to the observed scatter. Correcting for the con- 
tribution from the uncertainties implies an intrinsic dis- 
persion in relaxation times of 0.3 dex. 

In order to look for correlations of quasar variability 
properties on luminosity, redshift, black hole mass, and 
E ddington rati o, we used the line ar reg r ession method 
of [KelH (|2007[) . The method of iKelivl (pOOTi ) takes a 
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Fig. 5. — Distribution of the best fit CAR(l) process parameters 
for the 100 quasars in our sample. The characteristic time scales of 
AGN optical light curves are 10 < r < 10'* days, the amplitudes of 
short time scale variations are a < 0.02 mag day~*/^, and the am- 
plitudes of long time scale variations are < 40%. The uncertainties 
on the characteristic time scales are large, and the true dispersion 
in T is likely ~ 0.3 dex. 



Bayesian approach to linear regression, and accurately 
accounts for intrinsic scatter in the regression relation- 
ships, as well as measurement errors in both the depen- 
dent and independent variables. Measurement errors can 
be large for both the estimates of t and Mbh, and thus 
can h ave a significant effect o n the observed correlations 
(e.g.. lKellv fc Bechtoldl[2007l) . Therefore, it IS necessary 
to correct for the measurement errors when attempting 
to recover any underlying trends. 
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TABLE 2 

Results from CAR(1) Process Fits to Quasar Light Curves 



RA 


DEC 


logr^ 


Conf. Int. for logr 


log 


Err. in logcr 


J2000 


J2000 




95%, day 


mag day^/'^ 




00 29 13.7 


+13 16 03.9 


3.44 


[2.59,6.17] 


-1.97 


0.05 


00 47 15.8 


-72 41 12.2 


C 






00 49 34.4 


-72 13 09.0 


2.97 


[2.33,5.50] 


-1.95 


0.05 


00 51 16.9 


-72 16 51.1 


2.10 


[1.84,2.70] 


-1.60 


0.03 


00 54 52.1 


+25 25 39.0 


2.71 


[2.18,4.82] 


-1.81 


0.05 


00 55 34.2 


-72 28 30.0 


2.56 


[2.08,4.14] 


-2.31 


0.06 


00 55 59.6 


-72 52 45.1 


3.07 


[2.38,6.86] 


-1.85 


0.02 


01 01 27.8 


-72 46 14.4 


3.04 


[2.39,5.22] 


-2.04 


0.04 


01 02 14.4 


-73 16 26.8 


2.75 


[2.22,5.28] 


-2.23 


0.05 


01 02 34.7 


-72 54 22.2 


2.41 


[2.04,3.96] 


-1.80 


0.03 


01 07 21.7 


-72 48 45.8 


2.34 


[2.00,3.74] 


-1.89 


0.03 


04 46 11.0 


-72 05 09.0 


3.03 


[2.41,5.34] 


-2.20 


0.05 


04 53 56.6 


-69 40 36.0 


1.78 


[1.60,2.07] 


-1.93 


0.02 


04 56 14.2 


-67 39 10.8 


3.26 


[2.55,5.53] 


-2.17 


0.06 


05 00 17.6 


-69 32 16.3 


c 







Note. — Table [2] is published in its entirety in the electronic edition of the Astro- 
physical Journal. A portion is shown here for guidance regarding its form and content. 
^ The logarithm of the characteristic time scale of the quasar lightcurve, in days.'' The 
logarithm of the standard deviation in the input process to Equation Q. The standard 
deviation of quasar flux variations on time scales of 1 day is expected to be equal to 
a.'^ The CAR(l) process provided a poor fit to these lightcurves, and the data were 
not used in our regression analysis. 



4.1. Dependence of Quasar Variability on Luminosity 
and Redshift 

In order to investigate whether quasar variability prop- 
erties depend on luminosity, redshift, or both, we per- 
formed a regression analysis. Throughout this section, 
the luminosity will always be taken to be \Lx at 5100A. 
There has been considerable debate over whether quasar 
variability is correlated with luminosity or redshift, and 
the artificial correlation between the two has made it dif- 
ficult for previous work to uncover the true intrinsic cor- 
relation. However, because we perform a linear regresion 
of variability properties on both L and z simultaneously, 
we are able to break the degeneracy between L and z. 
This is because the multiple linear regression describes 
how variability depends on L at a given z, and likewise 
for z at a given L. 

In Figure[6]we show the relaxation time r as a function 
of luminosity and redshift, and in Figure [7] wc show a as 
a function of luminosity and redshift. The results of the 
regressions for r are 

logr = (-10.29 ± 3.76) + (0.29 ± 0.08) log XLx [da^) 
logr = (2.32 ± 0.10) + (1.12 ± 0.41) log(l -f z) [da](§!]l) 
logr = (-8.13 ±0.12) + (0.24 ±0.12) log ALa 

+ (0.34 ± 0.58) log(l + z) [days], (22) 

(23) 

and the results of the regression for a are 

log cr^ = (4.73 ± 2.34) - (0.19 ± 0.05) log ALa [R magV##] 

log cr2 = (-3.84 ±0.06) 

- (0.32 ± 0.25) log(l + z) [RmagVday] (25) 
log = (8.00 ± 3.29) - (0.27 ± 0.07) log XLx 

± (0.47 ± 0.33) log(l + z) [RmagVday]. (26) 

(27) 

There is a statistically significant correlation between r 
and both L and z, where the light curve relaxation time 



increases with increasing L and z. In addition, there is 
a statistically significant anti-correlation between a and 
L, implying that the short time scale variance decreases 
with increasing L. However, the multiple regression re- 
sults show that the luminosity trends are the dominant 
ones, and that there is no significant trend between either 
T or cr and z, at a given L. This therefore implies that 
the observed trends with redshift are caused by the artifi- 
cial correlation between L and z resulting from selection 
effects. 

We also looked for trends of the light curve variance, 
Ta'^/2, with L and z and found statistically significant 
evidence for a correlation between the variance of the 
light curve and z. Both the Spearman and Kendall rank 
correlation statistic was significant at 3cr, although there 
is considerable scatter in the correlation. This correlation 
is most likely a reflection of the well-known fact that 
quas ar emission at shorter wa velengths is more variable 
(e.g., IVanden Berk et al.ll2004D 

4.2. Dependence of Quasar Characteristic Time Scale 
and Variability on Mbh 

We also looked for trends in quasar variability prop- 
erties with Mbh and the Eddington ratio, L/LecIcI- In 
this work wc assume a constant bolo metric correction o f 
Cboi = 10 to the luminosity at 5100A (jKaspi et al.ll2000l ). 
However, we stress that a constant bolometric correc- 
tion can introduce significant error in the Eddington ra- 
tio (e.g., B^asudcvan & Fabian 2007; Kelly et al. 2008), 
and that our use of Cboi = 10 is only suggestive. Strictly 
speaking, what is being used in the following regressions 
is the ratio ALA(5100A)/Ms_f/, and in general this will 
not equal the true Eddington ratio. In Figure [S] we also 
show T as a function of Mbh and L/ Lndd, and in Figure 
[7] we also show cr as a function of Mbh and L/ Lndd- 

The results of the regressions for r are 

logT = (-2.29 ± 1.17) + (0.56 ± 0.14) logA/sH [da^) 
log r = (2.50 ± 0.24)- (0.06 ±0.27) log L/Lfeoi [daf^!^) 
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log Mujj / MgjjN Estimated log L / L^^^ 

Fig. 6. — The characteristic time scale of the optical light curves for the AGN in our sample as a function of optical luminosity, redshift, 
black hole mass, and estimated Eddington ratio. For clarity, we only show error bars for a random fraction of the data points in the top two 
panels, and we only show the error bars for the sources with Mgu estimated from reverberation mapping in the bottom two panels. The 
straight lines denote the best fit linear regression. There is a significant trend for t to increase with increasing Mgg, and less significant 
trends between t and \Lx or z. There is no significant trend between t and L/LEdd- 



and the results of the regression for a are 

logcr2 = (0.33 ±0.73) 

-(0.52 ±0.08) log Mb// [R magVday] (30) 

logcr2 = (-4.33 ±0.19) 

- (0.25 ± 0.22) log L/Lsdd [R magVday].(31) 
Based on these regression results, there is a statistically 
significant correlation between r and Mbh, where the 
relaxation time increases with increasing Mbh- In addi- 
tion, there is significant evidence that a decreases with 
increasing Mbh- However, there is no evidence for a 
dependence of r or cr on the Eddington ratio. 

In this work we ha ve employed the linear regression 
technique described in lKellvl (|2007l ). However, this tech- 
nique assumes symmetric errors in r, but in reality the er- 
ror distribution is asymmetric. In order to assess whether 
this has any ef f ects o n our results, we modifed the tech- 
nique of iKellvl ()2007[ ) to handle asymmetric errors in r, 
and used our modifed code to recompute the dependen- 
cies involving r. However, we did not notice any differ- 
ence in the slopes by assuming symmetric or asymmetric 
error bars. This is most likely because of the assump- 
tion of a linear relationship between log r and the other 
quantities. The asymmetric error bars for the high-mass 
quasars imply that the variability time scales for these 
sources are consistent with time scales much longer than 
the span of the lightcurve, which would lead to a steep- 
ening in the slope. However, this steepening of the slope 
would be inconsistent with the (better determined) time 



scales for the low-mass AGN, since a steeper line would 
predict shorter time scales than are observed. So, the 
assumption of a linear relationship, in combination with 
the well determined time scales at the low-mass end, 
counteracts the effects of the asymmetric error bars at 
the high-mass end. However, if we were to assume a non- 
linear form for the dependence of logr on log Mbh, such 
that the slope increases for the high mass objects, then 
the asymmetric error bars would likely have a noticeable 
effect, since the low-mass objects would no longer convey 
a lot of 'information' about the regression relationship at 
the high mass end. 

Due to the correlation between L and Mbh, it is un- 
clear whether the observed dependency of r and a on 
these quantities is real for both L and Mbh, or whether 
one correlation is simply a reflection of the other. Sim- 
ilar to breaking the L-z degeneracy, we can investigate 
which correlation is the fundamental one, or if both are, 
by performing a multiple regression of t and cr on L and 
Mbh- The results are 



loga^ = (-3.83 ± 0.17) - (0.09 ± 0.19) log 



10^5 erg s-i 



(0.25 ±0.24) log (^^^1^ ) [RmagVday] (32) 



(80.4^^6.9) 



/ Mbh 
1045 erg s-i 



-0.42±0.28 
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Fig. 7. — Same as Figure|6] but for the variance in the short time scale variations, . There is a significant trend for a to decrease with 
increasing Mbh, and a similar but less significant trend between a and XLx- 
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Here, we have expressed t as a function of L and Mbh in- 
stead of log T for more direct comparison with the charac- 
teristic time scales described by Equations (I3|)-(l5]). The 
joint probability distributions of the slopes are shown in 
Figure [5] for both regressions. It is unclear whether a 
depends on solely Mbh, solely L, or both Mbh and i, 
although the data favor a dependence on Mbh over one 
on L. However, there is significant evidence that the 
relaxation time scale depends on at least Mbh, and pos- 
sibly on L as well. In fact, the dependence of r on Mbh 
has steepened, and the relationship described by Equa- 
tion psp is similar to that for the orbital or thermal time 
scale, assuming a viscosity parameter of a ~ 10^'^ and 
location of the emitting region to be Rs ^ 100. 

5. DISCUSSION 

5.1. Comparison with Previous Work 

Most previous work on quasar optical variability has 
been based on analysis of structure functions or power 
spectra, either of individual quasars or an ensemble of 
quasars. From these studies, an anti-correlation be- 
tween varia bility and luminosit y has often emerged (e.g.. 
Hook e t aiJll994l:lGarcia et al.lll999HVanden Berk et all 
2004; WilhiteIirS][2QQl), while results on a variability- 
redshift correlation have remained mixed. Our result 
that long-term quasar variability is uncorrelated with lu- 
minosity may appear to be in conflict with previous work. 
However, the evidence for a variability-luminosity cor- 
relation is considerably weaker in the studies that have 



computed variability measures for individual objects. In- 
deed, studies that have compared variability with lumi- 
nosity for individual objects have noted the signiflcant 
scatter in the relationship, producing a very weak corre- 
lation, often leading to a detection of 'moderate' statis- 
tical significance at best. Ensemble studies, on the other 
hand, cannot investigate the scatter in the relationship 
and therefore cannot assess the strength of the correla- 
tion. Instead, ensemble studies simply look for an av- 
erage trend in variability properties, and, given enough 
quasars, are able to detect even a weak trend of vari- 
ability with luminosity. Furthermore, quasars exhibit a 
range in characteristic time scale and variability ampli- 
tude at a given luminosity or black hole mass, and it is 
unclear how this affects an ensemble structure function 
or power s pectrum. 

Recently IWold et all (|2007D have reported a correla- 
tion between optical variability and Mbh on time scales 
t > 100 days, but no correlation is seen on shorter time 
scales. In contrast, we observed an anti-correlation be- 
tween Mbh and short time scale variability, but no corre- 
lation betw e en M bh and long time scale variability. The 
IWold et al.l (|2007l) res ult is u nexpected, since Mbh and 
L are correlated ie.s... iPeter son et al. 2004), and previ- 
ous studies have found that variabilit y is anti-cor r elated 
with L, even on long time scales. The lWold et~aLl (|2007D 
result is based on an ensemble structure function, and 
the uncertainties on the structure function for the high 
Mbh bins are large. F urthermore, the time sampling 
of the IWold et aH (|2007f ) sample for the high M bh bins 
is worse than for the low Mbh bins, and windowing ef- 



AGN Variability 



11 




1.0 
0.5 



& 0.0 



t/2 

X 

m 



-0.5 

-1.0 
-1.5 



-1.0 



Regression 


- 







-0.5 



0.0 



0.5 



L slope 



Fig. 8. — Posterior probability distribution for the values of the coefficients in a linear regression of the characteristic time scale o f qu asar 
optical variations, r, (left) and the magnitude of short time scale variations, a, (right) as a function of Mgu and XLx (see Eg. 1321 and 
1331 ). The contours correspond to approximate 50%, 75%, and 95% joint confidence regions. While there is significant evidence that t 
depends on at least Mgu, it is unclear if there is an additional dependence on luminosity. In addition, while it is clear that a depends on 
either Mbh or \L\, it is unclear whether the dependency is on Mbh, ^L\, or both. 



fects due to the finite length of the time series may be 
at work here. When consid erin g v ariab ility measurement 
of individual sources, Wold et al,, (2007) still find a posi- 
tive correlation between variability and Mbh- However, 
while this correlation is statistically significant, it is very 
weak and exhibits considerable scatter. We performed 
a Kendall and Spearman rank correlation test between 
long time scale variability and the estimated black hole 
mass, and also find a marginally significant correlation, 
but the significance disappeared when we accounted for 
the uncertainty in the mass estimates. In addition, be- 
cause the long time scale variability for a CAR(l) process 
is Tcr^/2, Equations ([28|) and ([30| imply that the long 
time scale variability depends only weakly on Mbh , if at 
alL 

ICoUier fc Petersonl (|2001h calculated structure func- 
tions of optical light curves for 12 low-z AGN with re- 
verberation mapping data. Consistent with our work, 
they find a correlation between Mbh and characteris- 
tic time scale, where a characteristic time scale was de- 
fined as the location of a break in the structure fun ction. 
The time scales found bv lColher fc Petersonl (|2001h were 
consistent with dynamical or disk thermal time scales, 
although they tended to be somewhat shorter than those 
derived in this work. This difference may be explained 
by somewhat different definitions of 'characteristic' time 
scale between their work and ours, and systematic er- 
rors in the estimated structure functions caused by time 
sampling effects. 

Two of the sources in our sample, NG C 5548 and 
NGC 4151, w ere st udied by ICzernv et all (|1999[ ) and 
ICzernv et al.l ()2003[ ). respectively, using spectral tech- 
niques. These authors fit a functional form to the optical 
power spectra similar to that implied by the CAR(l) 
model, with the exception that they allow the high- 
frequency slope to vary. In contrast, the CAR(l) model 
we employ fixes this slope to be —2 i.e., Px{f) oc 1/ P- 
For NGC 5548, ICzernv et al] (fT999l) find evidence for a 
flattening of the optical power spectrum on time scales 
logT — 3.2 ± 0.3 days. This is longer than our observed 
characteristic time scale of logr = 2.3 ± 0.16 at w 2.5(7 
significance. This difference is only marginally signif- 



icant, and may be explained by the slightly different 
parameteric forms assume d for the p ower spectra, the 
longer light curve used bv ICzernv et al. (1999), the dif- 
ferent photometric bands, and errors introduced by the 
win dow function and sm oothing of the power spectra in 

the ICzerny et al.l (I1999D analysi s. 

For NGC 4151. ICze rnv et al.l pOOl did not find any 
evidence for a flattening to white noise of the _B-band 
optical power spectra over the time scales probed by their 
90 yr light curve. However, they flnd some evidence that 
the power spectra flattens below time scales of t ~ 100 
days, in agreement with the characteristic time scale we 
find of logr = 2.2±0.78 days. Furthermore, we note that 
our analysis constrains the characteristic time scale for 
NGC 4151 to be < 30 yr at 99% confidence. If the power 
spectra does flatten on time scales ~ 30 yr, it would be 
difficult to de tect in the spectral analysis employed by 
ICzernv et al.l ([2003), as such a turn-over would occur in 
the few lowest frequency bins. This is especially true 
when one considers that the lowest frequency bins of the 
power spectra are also among the bins most biased by 
windowing effects and smoothing. 

5.2. Connection with Accretion Physics 

In this work we have found that both the character- 
istic time scale of quasar light curves, and the magni- 
tude of the short time scale variations, depend on black 
hole mass. This, in combination with the evidence from 
reverberation mapping, strongly argues that the source 
of quasar optical variability is intrinsic to the accretion 
disk. In addition, we have found that the characteris- 
tic time scales of quasar light curves are similar to what 
would be expected for disk dynamical or thermal time 
scales, assuming a viscosity parameter of a ~ 10""^. Re- 
cent MHD simulations of radiation-dominated AGN ac- 
cretion disks have found that the thermal time scale is 
shorter than t hat implied by the standard a-prescription 
(|Turneill2004l ). making the association of r with thermal 
time scales more consistent. Our results imply that on 
time scales shorter than torb or tth, the accretion disk 
has difficulty generating variations in optical flux in re- 
sponse to random variations of some input process, such 
as, for example, a time varying magnetic field. Instead, 
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these short time scale variations get 'smoothed out', cre- 
ating a 1//^ power spectrum for frequencies higher than 
1/torb or ~ l/tth- 

While the quasar characteristic time scales are consis- 
tent with both orbital and thermal time scales, we find 
it more appropriate to associate these time scales with 
tfh, as we might expect some sort of periodic activity 
in the light curves if the flux variations were driven by 
orbital motion. In addition, quasar optical emission is 
thought to be thermal emission fro m an optically thick 
accre tion disk (e.g., Krolik 1999. : Frank. King, fc Raind 
Bool , and quasars tend to be bluer as th ey brighten 
(;e.g.. lGiveon et al.lll999l : lTrevese et al.ll200lD . suggesting 
that the flux variations are due to thermal variations. 
An association with orbital time scales cannot be ruled 
out, and indeed, the characteristic time scale of the tur- 
bulence driving the magnetic energ y and kinetic energ y 
density in the disk is ~ tort (e.g., iHirose et al.l 120081 ). 
However, if the radiation energy density in the disk is 
due to dissipation of magnetic and kinetic energy den- 
sity, then it is expected to respond to fluctuations in the 
magnetic and kinetic energy densities with a character- 
istic time scale of order the thermal time. 

In order to interpret the CAR(l) process in terms of 
accretion disk physics, we rewrite Equation ([T]) as 

d\ogL{t) = - -{log L{t)-^i)dt+(T[B{t+dt)-B{t)]. (34) 
r 

Here, L(t) denotes the luminosity of the quasar at time 
t, ij, = rb is the mean value of the quasar light curve, 
and B{t) denotes Brownian motion. In Equation ()34|) we 
have used the fact that the derivative of Brownian mo- 
tion is white noise, i.e., e{t) = dB{t) = B{t -\- dt) — B{t). 
Brownian motion is a non-stationary random walk pro- 
cess that has a power spectrum P{f) oc 1//^, and is 
described by Equation ([1]) in the limit r — s- oo. In ad- 
dition, Equation (I34p implicitly assumes that the vari- 
ance in the random variable dB{t) = B{t + dt) — B{t) 
is Var{dB{t)) — dt. The notation B{t) should not be 
confused with the common physics notation of denoting 
the amplitude of a magnetic field; although in this work 
we suggest that B{t) may be associated with a turbulent 
magnetic field, the term B{t) should be understood as 
referring to a Brownian motion. 

Writing the Equation for a CAR(l) process as Equa- 
tion (|34p reveals a number of interesting properties of 
this process. First, we note that the first term on the 
right side is what keeps the time series stationary. Con- 
sidering only this term (i.e., a = 0), d\ogL{t)/dt is neg- 
ative when the value of L{t) is brighter than the mean, 
and d\ogL{t) I dt is positive when L{t) is fainter than the 
mean. Therefore, the first term on the right side stabi- 
lizes the process by always driving L{t) toward its mean 
value, while the second term generates random perturba- 
tions to d\ogL{t)/dt that cause L{t) to deviate from its 
expected path. For highly accreting objects like quasars, 
the accretion disks are expected to be radiation pressure 
dominated in the regions emitting the optical and UV 
flux. Under the standard a-prescription for the viscosity, 
where the viscous torque is assumed to be proportional 
to the total pressure, a radiation pressure dominated disk 
is unstable to perturbations in the he ating rate (e.g., 
iShakura fc SunvaevI fl976t iKrolikI [l999l) . The fact that 
the quasar light curves in our sample are described well 



by a CAR(l) process with relaxation times similar to 
disk thermal time scales rules out instabilities in the 
disk that grow as ~ tth , consistent with res ults obtained 
from 3-dimensiona l MHD simulations (e.g.. lTurneHl200'l : 
'Hirosc et al.' 2008*). 

From Equation (p4|) it is apparent that the stochas- 
tic input into the differential equation, which drives the 
random variations in L, is itself a stochastic process. In 
particular, a random deviation in the input process, B{t)^ 
over a time interval dt causes a random perturbation to 
the change in logL(i) expected over the interval dt, with 
a controlling how sensitive d\ogL{t) is to dB{t). For 
example, the input process could be due to variations in 
accretion rate, or perturbations caused by a time- varying 
magnetic field. Indeed, recently some authors have mod- 
eled quasar X-ray variability as being driven by varia- 
tions in a magnetic field, with the magneti c field den- 
sity being modelled as an AR (1) process (e.g..lKing et al.l 
Ila04; Mayer & Pringlc 2006t I.Taniuk & Czernvll2007i r 

Hir ose et al.l (J2008) have argued for a similar in- 
terpretation of quasar optical variability, based on 3- 
dimensional MHD simulations. They argue that fluctu- 
ations in the magnetic field are dissipated, transferring 
energy from the magnetic field to heat in the plasma, 
creating thermal fluctuations in the accretion disk. The 
fluctuations in the heat content of the disk then create 
fluctuations in L(i); however, because the disk radiation 
cannot react to changes in heat content on time scales 
less than the thermal time scale, the shorter time scale 
fluctuations in flux are smoothed out and damped. Short 
time scale variations in radiation are correlated because 
the radiation energy density of the disk has not had time 
to completely react to the change in heat content induced 
by the change in magnetic energy density. However, on 
time scales t > tth^ the disk has had time to adjust to 
the heat content variations, thereby 'forgetting' about 
the previous perturbations in heat content. The result is 
a red noise power spectrum on time scales t <tth, and a 
white noise power spectrum on time scales t >tth- 

Within this interpretation, the parameter a describes 
the amplitude of variations in flux caused by variations 
in the magnetic energy density. Therefore, the anti- 
correlation between a and Mbh may imply two things. 
First, it may imply that the magnitude of the stochastic 
input process, B{t), possibly associated with a turbu- 
lent magnetic field, decreases with increasing Mbh- Or, 
it may imply that the sensitivity of changes in radia- 
tion energy density to changes in the turbulent magnetic 
field depends on Mbh, where fluctuations in the mag- 
netic fleld create smaller fluctuations in L(i) as Mbh 
increases. Alternatively, it could imply that both effects 
are at work. 

5.3. Microvariability and Implications for the 
Radio-loud/ Radio- quiet Dichotomy 

Modelling quasar variability as a stochastic process 
provides an opportunity to unify both short and long 
time scale variability as the result of a single process. 
The source of variations in radio-quiet quasar optical lu- 
minosity over time scales of hours (so-called 'microvari- 
ability' or 'intranight variability') has remained a puz- 
zle, although reprocessing of X-ra ys or a weak blazar 
component have been suggested (jCzernv et al.l [200& ). 
In general, microvariability in radio-quiet quasars is 
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not detected above t he photometric un certainty (e.g., 
iGupta fc Jo"shil l2005t ICarini et al.l l2007f ): however, for 
those sources for which it is detected the standard devia- 
tion i n the variabihty over the course of a night is ^ 0.01 
mag (iGopa l-Krishn a et ani2003t IStaUn et al.l l2004'. '2005|; 
iGupta fcJoshi 20051 ). As noted in the Appendix, for 
a CAR(l) process the standard deviation of short time 

scale variations is ~ aVKi. As can be seen from figure [H 
our best-fitting CAR(l) processes for quasar light curves 
predict variations of < 0.02 mag over 8 hours, consis- 
tent with what has been observed for radio-quiet quasars. 
Furthermore, it implies that the amplitude of microvari- 
ability should decrease with increasing black hole mass. 

To further investigate whether our best-fitting CAR(l) 
processes are consistent with observed microvariability 
of radio-quiet quasars, we calculate the expected distri- 
bution of microvariability amplitude for three different 
values of black hole mass: Mbh = 10^, 10®, and IO^Mq. 
We simulate light curves using the values of r and a cal- 
culated from Equations ([28]) and ([30| for each value of 
Mbh- Each simulated light curve had 40 data points, 
regularly sampled over the course of 5 hrs; these values 
were chosen to be con sistent with recent observat ions of 
micro variability (e.g., IStahn et al.l 120051 : [Gupta fc Joshil 
|2005() . For simplicity we ne glect observational error. Fol- 
lowing |Romero^r^ ()1999f ) we calculate the amplitude 
of microvariability from the difference in the maximum 
and minimum observed flux values, expressed as a per 
cent: 

Ip ^ == {Lmax{t) - Lrmn{t)) ■ (35) 



Here, L{t) is the average value of the light curve. 
The results are shown in Figure [9l As can be 
seen, the amplitudes of microvariability expected from 
our best fit CAR(l) processes are consistent with 
what is observed, at least f or radio-quiet quasars 
fe-g.-IGopal-Krishna et an[200a IStalin et al.l[200l 120051 : 
iGupta fc Joshill2005ri . 

Radio-loud quasars, on the other hand, are known to 
exhibit stronger microvariability, and the intranight fluc- 
tuations are often a bove the detection threshold (e.g., 
iGupta fc Joshi|[2005l ). If the physical mechanism for mi- 
crovariability is the same in both radio-quiet quasars and 
radio-loud quasars, then the observational result that 
radio-loud quasars exhibit stronger microvariability im- 
plies a correlation between a and radio- loudness. To 
test this, we compare the values of the radio-loudness, 
Recm, for the lGiveon et all (|1999D and AGN Watch sam- 
ples. The radio-loudness parameter, Recm, is defined 
in the standard way to be the ratio of flux density at 
6cm to the flux density at 4400A; the traditional divi- 
sion between radio-quiet and radio-lo ud sources occurs 
at Recm — 10 (|Kellermann et al.lll989D . Figure [TOl shows 
the best-fit values of cr^ , which give the optical variance 
on time scales of w 1 day, as a function of radio-loudness. 
As can be seen, a is anti-correlated with radio-loudness, 
showing that our best-fit CAR(l) process parameters 
predict microvariability should actually be reduced in 
radio- loud quasars. This is inconsistent with the obser- 
vational result that microvariability is stronger in radio- 
loud quasars, and implies an additional variability mech- 
anism in radio-loud objects on time scales shorter than 
those probed by our study. 



1.5 




Microvariability Amplitude, R-band [%] 

Fig. 9. — Expected probability distribution of radio-quiet quasar 
microvariability amplitude, ip, predicted from the best-fit CAR(l) 
process parameters, for Mbh = IO^Mq (solid line), Mbh = 
10* M0 (dashed line), and Mbh = 10^ (dashed-dotted line). 
The parameter ip is defined to be the maximum observed difference 
in flux values observed over a 5 hour period, in per cent. The am- 
plitude of microvariability predicted from our best-fitting CAR(l) 
processes is consistent with what has been observed in radio-quiet 
quasars, implying that the same process drives both the long and 
short time scale variations in these objects. 
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Fig. 10. — The best-flt values of ct as a function of radio-loudness, 
Recm- For the CAR(l) process assumed in this work, short-time 
scale variations have a dispersion of ~ a\/dt on time scales 
dt. The amplitude of microvariability decreases with increasing 
radio-loudness, opposite the trend found from direct observations 
of microvariability in radio-quiet and radio-loud sources. This, 
along with the fact that our study only probes time scales > 2 
days, implies that radio-loud objects have an additional variability 
component that operates on intranight time scales. 



Our analysis suggests that radio-loud quasars have an 
additional variability mechanism that operates on times 
scales < 1 day. Because the time sampling for the 
lightcurves analyzed in this sample is > 2 days, we are 
only able probe variations on time scales 2 days ^ St < 
7 yrs when fitting the CAR(l) process. Therefore varia- 
tions on time scales < 1 day are not used in estimating a. 
The fact that the magnitude of microvariability predicted 
by the best-fit CAR(l) processes for radio-quiet quasars 
is consistent with direct observations of microvariability 
implies that the microvariability in radio-quiet quasars is 
caused by the same process that drives the longer time- 
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scale variations. In this work we have suggested that this 
process, which is well-described as a CAR(l) stochastic 
process, is the result of a turbulent magnetic field driv- 
ing thermal fluctuations in the accretion disk. If this 
process is also the source of micro-variability in radio- 
loud objects, then the anti-correlation between a and 
radio-loudness would imply that microvariability should 
be reduced in radio-loud quasars. However, because the 
opposite is true, this implies that radio-loud objects have 
an additional source of variability on time scales < 1 day, 
but does not operate on the time scales probed by are 
study, i.e., 2 days 5t <1 yrs. This additional source 
of optical variability in radio-loud objects is likely due 
to the presence of a relativistic jet, as has often been 
suggested. 

Because the magnitude of microvariability predicted by 
our best-fit CAR(l) process is consistent with observa- 
tions of microvariability in radio-quiet quasars, it is not 
necessary to invoke an additional physical mechanism to 
explain the short time scale variations in these objects. 
Instead, the CAR(l) process is able to explain both short 
and long time scale variations as being driven by the 
same process, possibly thermal fluctuations being dr iven 
by a turbulent magnetic field. ICzernv et al.l ()2008D in- 
vestigated three different models for microvariability in 
radio-quiet quasars and concluded that the most promis- 
ing source of the m icrovariability was a weak blazar-like 
jet. Furthermore, ICzernv et atl (|2008D 



concluded that 



luminosity fluctuations driven by the magneto-rotational 
instability are not strong enough to create microvariabil- 
ity. Our e mpirical result s are in consistent with the con- 
clusions of ICzernv et al.l (I2008D in that we do not find 



any evidence for an additional process that drives the 
microvariability. This is not to say that our results im- 
ply that a weak jet does not exist in radio-quiet quasars, 
but rather that if a weak jet does exist then it's variabil- 
ity amplitude is small compared to that caused by the 
process that also drives the longer time scale variations, 
such as a turbulent magnetic field. 

It is u nclear how sensitive the conclusions of 
ICzernv eT al. (2008) are to the specific s of th eir model, 
and the investigation of ICzernv et al.l ([2008) does not 
necessarily rule a turbulent magnetic field as being the 
source of microvariability in radio-quiet quasars. Of par- 
ticular interest is the manner in which th e time depen- 
dence of the mag netic field is modele d. ICzernv et al. 
2008fl . as well as EEg et a l. (2004), IMaver k. Pringk 
2006( ). and iJaniuk fc Czernv (2007.), have modeled the 



time series of the magnetic field as an AR(1) process: 



Ui+i = aARUi + Ci 



(36) 



where the amplitude of the local magnetic field at the 
i*'' time step is proportional to Ui, and is a random 
variable uniformly distributed between ±-\/3. Because 
the characteristic time s cale of the magneti c field devel- 
opment is Tb kdTorb, (jCzernv et al.ll2008D used a time 
step equal to the orbital time. Ho wever, Equation ffl . 
in combination with Figure 10 of ((Czernv et al.ll2008f ). 
implies that the orbital time scale should be ~ 100-1000 
days, depending on the mass of the black hole. This re- 
sults in a time sampling of the magnetic field that is much 
longer than the time scales over which microvariability 
are measured, therefore artificially suppressing the am- 
plitude of microvariability in the simulated light curves. 



Moreover, the amplitude of microvariability can also be 
increased by increasing the variance of the random vari- 
able Ci. 

Following [King et all ([2004D . ICzernv et all ([2008[ ) as- 
sumed a value of aAR = 0.5. As described in the ap- 
pendix, a value of uar = 0.5 corresponds to a charac- 
teristic time scale for the equivalent continuous process 
of r = lAATorb, i-e., a value of kd = 1.44. A more ac- 
curate simulated light curve can be obtained by using 
a sampling interval in Equation (|36p that is small com- 
pared to that probed by microvariability, such as At ~ 
minutes. Then, the value of a can be related to the char- 
acteristic time scale of the magnetic field development as 
<^AR = s^p{—^t/ (kdTorb)}- In this case, the only free 
parameters are kd and the variance of e;. Further im- 
provement can be obtained by employing 3-dimensional 
MHD simulations to test whether a turbulent magnetic 
field can create microvariability with fluctuations ~ 2%. 

6. SUMMARY 

In this work we have modeled quasar light curves as a 
type of stochastic process called a first-order continuous 
autoregressive process. This statistical model has three 
free parameters: the characteristic time scale for the pro- 
cess to 'forget' about itself, r, the magnitude of the small 
time scale variations, a, and the mean of the time series, 
/i. We used this model to fit 100 quasar light curves 
at z < 2.8, including 70 quasars with black hole mass 
estimates. Our conclusions are summarized as follows: 

• Quasar optical light curves are often well described 
by a continuous autoregressive process (CAR(l)). 
For the quasars in our sample, the amplitude of 
variations on time scales long compare to their 
characteristic time scales (i.e., sim decades) are 
typically ^ 3-30%, consistent with previous work. 
The characteristic time scales of the quasar light 
curves vary between ~ 10 days and ~ 10 yrs and 
are consistent with accretion disk orbital or ther- 
mal time scales, assuming a viscosity parameter of 
a ^ 10""^ and location of the emitting region to 
be Rs ^ 100. In addition, the short time scale 
variations are < 0.03 mag day^^^^. 

• The characteristic time scales of quasar optical 
light curves are correlated with Mbh and lumi- 
nosity, while the magnitude of the short time scale 
variations are anti-correlated with Mbh and lumi- 
nosity. We did not find any evidence for an ad- 
ditional redshift correlation. A multiple regression 
analysis suggested that the primary correlation is 
with Mbh- At a given luminosity, the characteris- 
tic time scales depend on Mbh as 

-.J \ -0.42±0.28 
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1.03±0.38 



[days] 



where the errors are quoted at 68% confidence (Icr). 

• Within the CAR(l) process model of quasar light 
curves, the random perturbations to d\ogL{t) are 
caused by a stochastic process. This stochastic in- 
put process may be a turbulent magnetic field that 
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creates fluctuations in the accretion disk radiation 
energy density. Variations in the input process over 
an interval dt create variations in L{t) which are 
smoothed out on time scales shorter than an or- 
bital or thermal time scale. 

• The fact that radio-quiet quasar optical light curves 
can be well fit by a CAR(l) model suggests that 
it is not necessary to invoke an additional physi- 
cal mechanism to describe short time scale varia- 
tions in these objects. Instead, within the CAR(l) 
model, both short and long time scale variations 
are driven by an underlying stochastic process that 
causes random perturbations to d\ogL{t)/dt. This 
is supported by the fact that the intranight varia- 
tions predicted by our best fit CAR(l) processes are 
2% in the _R-band, consistent with observations 
of intranight variability. Furthermore, our best- 
fit CAR(l) parameters predict an anti-correlation 
between radio-loudness and microvariability, incon- 
sistent with what has been observed. This suggests 
that radio-loud quasars have an additional optical 
variability mechanism that operates on time scales 
shorter than those probed by our study, < 2 days. 

Before concluding, we stress that the CAR(l) model 
is a statistical model and not a physical model. 
Quasar light curves are stochastic in nature (e.g., 
ICzernv &: Lehtolll997| ). and the dependence of luminos- 
ity on the time- varying properties of the disk is complex. 
In this sense, the randomness in the stochastic model is 
not due to that fact that the physical processes them- 
selves are not deterministic, but rather is a reflection 
of our lack of knowledge of the complex physical pro- 
cesses that generate variations in flux. While a phys- 
ical model is needed in order to interpret the stochas- 
tic model in terms of accretion disk physics, and thus 
lead to a proper understanding of quasar light curves, 
the stochastic model is sufficient for modeling the data, 
given our current knowledge. Furthermore, much of the 
mathematical formalism of accretion physics is in the lan- 
guage of differential equations, suggesting that stochastic 
differential equations are a natural choice for modeling 
quasar light curves. 



The field of stochastic processes is a rich field with well- 
developed methodology, predominantly because of its im- 
portance in financial and economic modeling. We have 
utilized the CAR(l) model because of its simplicity, and 
because it allows us to perform statistical inference with- 
out having our results biased by the irregular sampling, 
measurement errors, and finite span of the time series. 
However, the CAR(l) model is the simplest of stationary 
continuous autoregressive processes, and additional flexi- 
bility may be achieved through the addition of higher or- 
der derivatives to Equation ([T]) . This provides a rich and 
flexible method of modeling the power spectra of quasar 
light curves without suffering from the windowing effects 
that can bias traditional Fourier and structure function 
techniques. For example, quasi-periodic oscillations can 
be modeled through the addition of second order deriva- 
tives to Equation (jT]) , as has been done in the analysis o f 
the frequency of sun spot numbers (jPhadke fc WulH974f) . 
In addition. Equation ([T]) can be generalized to a vector 
form, allowing the simultaneous modeling of quasar light 
curves across multiple observing wavelengths, and thus 
introducing additional constraints on physical models of 
quasar variability. Both the use of higher order terms and 
multiwavelength modelling of quasar light curves will be 
the subject of future research. 

A computer routine for fitting the CAR(l) model to a 
time series is available on request from B. Kelly. 
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APPENDIX 

DESCRIPTION OF AUTOREGRESSIVE PROCESSES 

The first order autoregressive process (AR(1)) process is a well-studied stochastic process (e.g., see IScarglel[l98lL 
and references therein), generated according to 

Xi = aARXi-i + (Al) 

where is a normally-distributed random variable with zero mean and variance cr^^, and the data Xi are observed 
at regular time intervals. The parameters of the AR(1) model are uar and cr^^, and uar is usually constrained as 
|c«a_r| < 1 in order to ensure stationarity; a time series is said to be stationary when its mean and covariance do not 
vary with time^. The case uar — 1 corresponds to a random walk. For quasar light curves, the n observed data points 
Xi, . . . ^Xn correspond to the observed fiuxes at times ii, . . . , t„, for ti = ti -\- {i — 1) Ai. 

The discrete AR(1) process is only defined for regularly sampled time series. However, astronomical time series 
are rarely regularly sampled, and often large gaps in time can exist. Furthermore, the AR(1) process is a discrete 
process, but the underlying physical process that gives rise to the observed fiux is continuous. Because of these two 
considerations, we instead model the quasar light curves as a first order continuous autoregressive (CAR(l)) process. As 

^ Technically, this is referred to as 'weakly stationary', while stationarity in the strictest sense implies that the probability distribution 
of the time series does not change with time. However, because we are only concerned with Gaussian processes in this work, the two 
definitions are equivalent 
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noted in §|3l the CAR(l) process is described by the following stochastic differential equation fe.g.. lBrockwell fc David 
[20nl) : 

dX{t) = --X{t)dt + aVdie{t) + bdt, T,a,t>0, (A2) 

T 

where, r is called the 'relaxation time' of the process X{t), and e(i) is a white noise process with zero mean and 
variance equal to one. Throughout this work we have assumed that the white noise process is also Gaussian. In the 
physics literature, Equation (|A2p is often referr ed to as an Orns tein-Uhlenbeck (0-U) process, and plays a central role 
in the mathemematics of Brownian motion; see IGillespid ([l99l for a review of the 0-U process. 
The solution to Equation (jA2p is 



A:(t) = e-*/^A:(0) + 6T(l-e-*/^) + cr / e-(*-^)/^dS(s). (A3) 



Here, ^(0) is a random variable describing the initial value of the time series, dB(s) is a temporally uncorrelated 
normally distributed random variable with zero mean and variance dt. Strictly speaking, dB{t) is an interval of 
Brownian motion, and the integral on the right side represents the stochastic component of the time series. From 
Equation (jASP it can be seen that if cr = 0, i.e., if there is no stochastic component, and if X{G) represents a random 
perturbation, X{t) relaxes to it's mean value with an e-folding time scale t; hence the identification of r as the 
relaxation time. If cr > 0, then the path that X{t) takes will vary randomly about the expected exponential relaxation. 
The expected value of X{t) given X{s) for s < i is 

E{X{t)\X{s)) = e-^'/^X{s) + 5t(1 - e"^*/^) (A4) 
and the variance in X{t) given X{s) is 

Var{X{t)\X{s)) = ^ [l - e'^^'/^] (A5) 

where At = t — s. If 2At/T ^ 1, then Equation (jASp implies that the variance on time scales much shorter than the 
relaxation time is « a^At. Therefore, can be interpreted as representing the variance in the light curve on short 
time scales. In addition, one can show that when the time sampling is regular with At = 1, then the CAR(l) process 
reduces to an AR(1) process with uar. — e^^/"^ and cr^^ = Ta^{l — e^^/'^)/2. 

In astronomical time series analysis it is common to interpret a light curve in terms of its autocorrelation function 
and power spectrum. The autocovariance function at time t' is defined to be the expected value of the product of X{t) 
and X{t + 1'), and the autocorrelation function is calculated by dividing the autocovariance function by the variance 
of the time series. The autocorrelation function of the CAR(l) process is 

ACF{t') = e-*'/^ (A6) 

Equation (jA6|l states that the correlations in C AR(l) light curv e fall off exponentially with lag t' , with an e-folding 
time equal to the relaxation time, r. Following IGillespid (|1995f) . the power spectrum of a process is computed from 
the autocovariance function rx(t') as 

/>oo 

Pxif) = 4 / rx{t') cos{2n ft') dt' , / > 0. (A7) 



For a CAR(l) process, rx{t') = [ra"^ /2)e * from which it follows that the power spectrum of a CAR(l) process is 

2ctV2 

The power spectrum of the CAR(l) process falls off as 1//^ on time scales short compared to the relaxation time, and 
flattens to white noise at time scales long compared to the relaxation time. 
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